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小 天 体 离 艇 元 仿真 软件 的 搭建 与 测试 
刘 沐 林 1 RBZ 
(南京 大 学 天 文 与 空间 科学 学 院 南京 210023) 

摘要 小 天 体 探测 是 当今 太阳 系 探测 的 一 大 热点 , 对 小 行星 演化 的 研究 有 助 于 人 们 了 解 太阳 系 的 起 源 . 演化 

研究 的 一 项 重要 内 容 是 小 天 体 结构 的 演化 , 即 小 天 体 在 多 种 力学 机 制作 用 下 自身 形状 与 结构 的 演化 . 在 小 天 

体 碎 石 堆 结构 的 假设 之 下 , 一 种 比较 常见 的 模拟 小 天 体 结构 演化 的 方法 是 离散 元 仿真 算法 (Discrete Element 

Method, DEM), 目前 国内 外 有 数 个 团队 开发 了 相关 算法 软件 . 本 文 介绍 了 团队 开发 的 《基于 DEM 仿 真 算法 的 

多 粒子 系统 模拟 软件 》 的 基础 理论 、 实 现 方法 以 及 加 速算 法 , 并 使 用 二 体 接 触 模型 、 声 速 传播 仿真 、 小 行星 内 

部 压强 、 小 行星 旋转 稳定 性 仿真 算 例 验证 了 算法 的 可 靠 性 . 
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中 图 分 类 号 : P132; ”文献 标识 码 : A 
1 背景 介绍 ee a 
HII E 点 H 是 直径 在 102- 104 mz. 间 的 小 天 体 H ae 门类 : Distinct Element Method, Modal ere 
是 碎 石 堆 结构 . 碎 石 堆 结构 最 早 被 用 于 地 质 学 , 如 Discontinuous Deformation Analysis, Momentum 
矿物 开采 叫 、 月 壤 研 究 四 等 随 着 人 们 对 于 小 行星 Exchange Method, 这 4 种 门类 对 应 了 不 同 的 接触 
性 质 认识 的 不 断 加 深 , 相关 研究 在 小 行星 领域 得 到 模型 和 算法 , 可 以 按照 接触 过 程 的 计算 方法 和 接触 
拓展 , 譬如 小 行星 自转 速度 上 限 回 、 小 行星 及 小 行 ”” 体 的 “ 软 硬 " 来 进行 划分 . 
星 族 的 磁 撞 起 源 &L6、 赶 星 结构 四 等 . 这 些 研究 的 Distinct Element Method (DEM) 即 本 文 工 作 

前 提 是 , 在 微 引力 环境 下 , 小 行星 的 组 成 是 在 自身 。” 所 使 用 的 方法 , 也 称 为 离散 元 仿真 算法 , 计算 中 

See ee a 忽略 接触 物体 本 身 发 生 的 形变 , 而 使 用 接触 体 之 
的 碎 石 堆 的 结构 , 而 基于 这 一 假设 所 得 的 研究 结果 间 的 相对 位 置 和 重 受 来 仿真 形变 的 过 程 ,， 即 硬 
目前 可 有 效 地 解释 地 面 和 就 位 探测 任务 的 一 些 观 。“ 球 的 软 碰撞 . 现 有 的 仿真 项 目 主要 有 : Parallel k-D 
测 现象 . Tree Gravity Solver (PKDGRAV)!"4J, Discrete El- 

在 小 天 体 碎 石 堆 结 构 的 前 提 下 , 可 以 使 用 离 。 ement Model Body (DEMBody) 21%, 


散 元 方法 来 对 小 行星 的 结构 与 演化 进行 


有 的 离散 元 仿真 方法 最 早 由 Cundalll 9 于 1971 年 提 
E, 用 于 解决 宕 体力 学 的 问题 . 随后 经 过 数 十 年 
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Modal Method 与 Distinct Element Method 相 
车 加 方法 来 模拟 接触 体 
F 紧密 堆积 的 物体 , 这 


方法 


65 4 


中 的 特征 模式 不 能 很 好 组 
方法 更 适合 在 下 松 的 系 
目 主 要 为 CICED3]. 


织 成 接触 的 边界 , 因此 这 
统 中 使 用 . 现 有 的 仿真 项 


Discontinuous Deformation Analysis, 即 著 名 


KW ANIES 
即 仅 计 入 
E 


接触 时 发 生 的 动量 
其 体 过 程 , 但 是 在 转化 过 程 中 将 接触 


卖 变形 分 析 方 法 . 在 计算 时 使 用 刚性 接触 ， 


LA 


-HH J 


进行 处 理 ， 


刚体 运动 和 弹性 
前 岩 体 力学 领域 被 大 量 使 用 的 计 入 
Momentum ees Method, 即 考虑 刚 


体 


本 作为 区 


转换 而 不 仿真 接触 的 


性 体 
的 形变 , 是 目 


方法 . 


体 的 


刚性 碰撞 , 接触 
分 离 继续 仿真 各 自 
开销 最 小 , 在 早期 的 离散 


为 开展 小 天 体 结构 演化 的 相关 研 
发 了 仿真 软件 《 
算法 的 多 粒子 系统 模拟 软 人 


离散 元 仿真 算法 开 


介绍 了 软件 构建 的 


便利 性 . 第 2 节 中 了 


介绍 了 仿 


的 运动 状态 . 这 一 方法 在 计 和 


上 


元 仿真 中 发 挥 重要 的 作 


真 软件 的 构 


基本 原理 、 使 用 的 


re 
学 模型 参数 的 选取 等 ; 第 3 节 介 


绍 了 使 


Paces 


并 对 结果 进行 了 讨论 


的 工作 进行 总 i 结 . 
2 ”基础 理论 


仿真 软件 


的 动力 学 模 


型 | 


F) (AstroDEM} 
节 以 及 测试 算 例 的 结果 , 经 过 
结果 比较 和 分 析 证 明了 方法 的 有 效 性 、 可 靠 1 


仿真 基本 单元 的 形 


R, 作者 基于 


基于 DEM 仿 真 


, 本文 


生 和 
ean, 包括 


力学 模型 和 力 


用 该 软件 的 
; 第 4 节 对 文章 


影响 , 不 同 的 颗粒 身 


准 静 止 角 等 性 质 . 
2.2 ” 非 球形 颗粒 单元 


虽然 球 


颗粒 单元 在 仿真 


元 形状 将 显 


学 模型 上 容 
o 


易 、 


境 下 , 颗粒 的 形 ; 


行 建 模 和 计 


J 能 会 导 
可 能 会 导 


E 


模型 [ 
BRR 

通过 不 同 

球 复 颗粒 


臻 整体 宏观 特性 
系统 


AR, FA 
筑 自 身 的 接触 力学 ， 
的 优点 . 对 于 复杂 的 非 球形 颗粒 结 


用 数 十 


构造 , 在 仿真 中 将 会 大 大 加 大 计 


的 构造 过 程 


中 ， 


E 


HI 


Nt 
$ 


影响 颗粒 堆 


中 


也 使 用 非 球形 颗粒 
真 . 非 球形 颗粒 形状 和 接触 模型 的 构 关 
4 和 多 面体 模型 上 5 9 两 种 方法 . 
pe Be AY DA BRIE UA A 
的 球形 颗粒 单元 县 加 


HER 
A 


aN 


一 步 需要 使 


LATE 
的 优点 , 但 
\ 上 共有 很 大 的 不 规则 性 ， 
的 显著 差异 ， 


AKA JI 
是 在 真实 


因此 , 一 


Na 


没 


元 来 进行 仿 


ERARIK 


元 作为 基本 单元 结构 ， 
的 方式 组 成 不 同 的 
颗粒 单元 的 接触 力学 来 构 
有 构造 便利 、 


i= 


显 奇 点 


有 日 


E, 


往往 需要 使 


至 数 百 个 球形 颗粒 单元 


EJ 


基本 单元 来 


H 


得 到 一 个 点 与 面 构 成 的 多 面体 形 ; 


同 算法 


[14, 17-19 


KEE, 再 使 用 球 艇 来 模拟 复杂 形 : 


和 质量 分 布 ; 第 三 步 使 用 数值 方法 计算 得 
型 的 惯量 张 量 , ; 
心 运动 进行 仿真 . 如 图 
在 不 同 精度 下 的 不 同形 ; 


造 的 球 艇 模 
转 、 扭 转 、 


We Bea AE 
真 模型 ， 


状 模 型 和 受 力 模型 组 成 , 这 两 者 目前 有 诸多 实用 模 
并 明确 本 文 仿真 软件 选 


型 可 参考 , 下 文 将 简单 介 台 


取 的 模型 . 
2.1 ”球形 颗粒 单元 


在 仿真 中 , 球形 颗粒 和 


元 是 最 简 上 


单 直 观 的 仿真 


i 


基本 单元 , 即 仿真 颗粒 为 具有 半径 7、 


Re 
wh 
ti 


度 p 的 匀 质 


球体 . 球体 颗粒 单元 不 仅 可 以 用 于 仿真 具体 的 仿真 


JE 


ce 


ae, 也 可 以 作为 下 面 球 簇 模型 的 基本 组 成 单元 ， 


模拟 具有 一 定 空间 复杂 度 的 非 球形 颗粒 形状 . 仿真 


颗粒 单元 的 形状 对 于 仿真 系统 的 动力 学 


生 质 有 较 


=. 


计算 机 软件 著作 权 登 


记号 : 2024SR0001744 


30-2 


多 


质心 


型 构造 的 


量 . 在 球 徐 模型 
i, 雷达 等 方式 


大; 第 


利用 点 和 面 在 区 域内 添 


KAY 


二 步 使 用 不 
加 球形 颗粒 
外 部 接触 面 


1[ 


到 所 构 


步 对 颗粒 碰撞 的 旋 


4 为 使 用 球 


看 体 模型 与 球 模型 


接触 力学 模型 


触 3 


天 颗粒 的 仿 


4 可 以 划分 为 面 接触 、 
形式 09 , 3 种 形式 的 接触 力 模型 均 作 


td 
Ji 


EH, 通过 其 质量 


BC HY ZH ae = fA h 
行 求 和 , 即 可 得 到 多 面体 旧 
以 及 惯量 张 量 民 


TE, XA h 
单元 所 受 的 力 与 力矩 
和 可 对 系统 进行 外 


体 所 


4 有 显著 的 不 同 . 一 方面 ， 
边 接触 、 角 接 


j 在 多 面 
有 面 元 的 受 


HE; 男 一 方面 , 搜索 入 


元 的 形状 的 复杂 性 


用 球形 颗粒 单元 所 使 用 的 仅仅 考虑 
离 的 搜索 算法 , 需要 考虑 


的 参数 . 


法 上 pal, | 


于 多 面体 颗粒 单 
E, 多 面体 的 接触 力学 不 能 完全 使 


两 颗粒 质心 时 
额外 的 多 面体 形状 等 更 多 


距 
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图 1 不 同 精度 参数 下 的 正八 面体 的 球 簇 模型 ( 左 、 中 ) 和 真实 颗粒 的 球 
RRDA (E). 左 列 展示 了 调整 球 簇 中 不 同 球体 之 间距 离 的 最 小 
值 dis 由 大 到 小 的 结果 , 此 时 随 着 din 的 减 小 , 正八 面体 的 表面 变 得 
光滑 ; 中 列 展示 了 调整 球 篮 中 最 小 组 成 球体 半径 rmin 由 大 到 小 的 结果 ， 
此 时 随 着 rmin 的 减 小 , 正八 面体 的 表面 细节 变 得 更 加 凸显 ; 右 列 展示 了 
两 个 参数 与 使 用 的 三 维 节点 数 最 大 比例 pmax 共 同 作用 在 一 个 真实 颗粒 
上 的 结果 , 在 不 同 的 参数 选择 下 ,颗粒 表面 的 精细 度 、 粗 糙 程 度 和 颗粒 

密度 分 布 有 所 不 同 D4 . 


Fig.1 Models of a regular octahedron (left, middle) and 
models of a real particle (right) under different spheres 
cluster precision setups. The left column displays different 
models of a regular octahedron with a decreasing value of 
the minimum distance dmin between each sphere in the 
cluster from the top to the bottom and it turns out that the 
surfaces of the models are becoming smoother; the middle 
column displays different models of a regular octahedron 
with a decreasing value of the minimum radius rmin of all 
spheres in the cluster from the top to the bottom and it 
turns out that the surfaces of the models are becoming more 
precise; the right column displays different models of a real 
particle with different din, Tmin, the maximum percentage 
of nodes used denoted as pmax and it turns out that every 
different setup shows different smoothness, precision and 


density distribution“). 


2.3 ”接触 力 模型 
接触 力 模 型 是 离散 元 仿真 的 核心 , 可 以 依据 颗 
粒 单 元 的 可 形变 性 划分 为 硬 球 和 软 球 模型 , 同时 ， 


也 可 以 依据 碰撞 过 程 的 数值 仿真 方式 划分 为 硬 接 
触 和 软 接触 . 这 两 种 划分 的 组 合 便 形 成 了 4 种 可 供 
使 用 的 离散 元 接触 力 模型 , 它们 具有 不 同 的 特点 ， 
可 以 适用 于 不 同 的 仿真 条 件 . 

在 本 工作 中 , 目前 主要 采用 球形 颗粒 单元 进 
行 计算 . 考虑 两 颗粒 间 的 接触 力 , 设 两 颗粒 1、2 的 
MAA AAT Me, TAMA ie RAT, Mv, 
AeA TAG, Ma, 颗粒 半径 分 别 为 Fr 和 rz， 
于 是 颗粒 1 受到 颗粒 2 的 接触 力作 用 和 力 抢 可 以 分 
别 写 成 : 


F =F, + Fr, 
M = -hù x Fr, 
其 中 , A 二 B® LARL, 接触 力 可 


|#—#3]” 
DAs AE LF AY) a Fi PS, 记号 n 代 表 法 向 
(Normal), 人 r 代 表 真 实 的 摩擦 力 (Real Friction). 我 
们 在 下 面子 节 中 详细 介绍 法 向 和 切 向 的 接触 力 模 
2.3.1 法 向 接触 模型 
考虑 两 颗粒 发 生 接 触 的 情形 , 对 于 非 刚体 颗粒 ， 
接触 时 将 会 发 生 形变 并 产生 接触 弹力 , 从 而 将 两 颗 
粒 向 相反 方向 弹 开 . 在 本 文采 用 的 方法 中 , 系统 不 
对 颗粒 的 形变 进行 仿真 , 而 是 通过 考虑 两 颗粒 之 间 
重 登 部 分 的 大 小 来 估计 接触 时 因 形 变 产生 的 弹力 
KANG, 得 到 两 颗粒 相互 重 又 部 分 的 大 小 后 , 即 采 
用 相应 法 向 接触 模型 和 切 向 接触 模型 对 接触 力 进 
行 仿真 . 图 2 为 两 个 接触 颗粒 接触 过 程 计 算 各 个 参 
数 的 示意 图 , 定义 两 颗粒 之 间距 离 为 d = |z 一 zl, 
1 和 /2 分 别 为 颗粒 1 和 颗粒 2 的 转动 臂 , 上 且 有 : 


Ne 2d 
pa anit 
2 2d , 


是 两 球体 的 重 登 量 , 用 于 表征 球体 的 形变 : 


E = ri +r — |£ 22|. 


1s $B i RO FL PL REA, 弹簧 模 
型 可 以 分 为 两 种 : RENEI REREN, 考 


65 卷 天 


虑 弹簧 的 弹性 项 和 阻尼 项 , 可 以 将 法 向 分 量 瓜 进 行 
如 下 分 解 : 


F, _ F, T Fy ) 
其 中 瓦 为 弹性 项 , 使 用 线性 弹簧 模型 可 以 表示 为 : 
F = kñ. 


Grain 1 Grain 2 


q2 ”两 颗粒 的 接触 过 程 计算 各 参数 示意 图 . 其 中 , z Za A 
中 心 之 间 的 距离 , 0: 和 12 分 别 为 颗粒 1 和 颗粒 2 的 转动 臂 ，6 是 两 颗粒 的 


Fig.2 Demonstration of parameters used in the simulation 
of two colliding grains. |%, — 72| is the distance between the 
two grains’ centers, lı and lz are rotating arms of grain 1 
and grain 2 respectively, and € is the overlapping distance of 


the two grains. 


AR ZG PE SSE BM TY] URRY: 
F, = kué” ñ , 
ku 为 球体 的 线性 弹性 系数 , kn 为 球体 的 非 线性 弹性 


A 
户 为 阻尼 项 , 线性 弹 得 


中 可 以 写作 : 


Fi = 一 yn ， 
在 非 线 形 弹 簧 模型 中 可 以 使 用 线性 弹簧 中 相同 的 
表示 , 也 可 以 写 为 B24; 
Fi = 一 ?55 th, 
7 为 阻尼 系数 ， 表征 球体 形变 的 速率 , A: € = Tip - 
Ah, HAG. SE: te = 0, — U+ x pi — We X po, 
Tp, = =r Də = Tof. 


2.3.2 ” 切 向 接触 模型 
XEHDET, 在 颗粒 表面 切 向 颗粒 受到 的 
弹力 可 以 表示 为 : 


Ë = —k,Dut , 


kE, 为 切 向 弹性 系数 , RA, = Fko, DE 为 切 
KEPLER, 定义 为 Dn = fy le ldt, 此 处 为 
次 接触 进行 总 时 间 , 读 为 切 向 速度 , 即 : 


I 


= 


U = tia — (Vig -A/)A, 


t 为 切 向 单位 向 量 , 于 是 有 t = (oh. 一 般 碰撞 时 会 取 
较 小 的 时 间 步 长 以 保证 仿真 的 准确 性 , 故 可 以 将 积 
分 表示 成 如 下 离散 化 的 形式 : 


Dy = 》 | 雹 | 上 

当 切 向 弹性 力 大 于 最 大 静摩擦 时 , 两 球体 之 间 
会 发 生 滑 动 , 此 时 切 向 力 表现 为 摩擦 力 , A: F = 
-ul Fli, 于 是 切 向 力 瑟 取 二 者 之 间 的 较 小 值 , 有 : 


Fi, = — min (ke Dy, uF Dê. 


于 是 两 颗粒 之 间接 触 过 程 中 施加 在 颗粒 1 上 的 
tea EM = fh x Fy. 

考虑 仿真 颗粒 与 固定 平面 之 间 的 作用 , 颗粒 质 
心 受 力 玉 和 颗粒 所 受 力 矩 邮 分 别 为 : 


F= F, + Fe, 
M = -hù x Fr. 


但 是 此 时 力 的 表达 式 有 所 不 同 . 首先 考虑 法 
aed pti a 

= yi, 但 是 此 时 有 & = ri- d, 其 中 必 为 
的 设 平面 的 标准 方程 为 Ax 
| + Oz 十 DD = 0, 颗粒 的 质心 坐标 (zx1, y1, 21); 
A A l, MEN RIAA: € 
to À = (Ù — rı xh) A= T- À, XFU E 
力 , 同样 有 : F = -k Dyf, (AEM EM RIAA 
所 不 同 . 考虑 切 向 速度 有 : 该 = 访 2 一 (tig -A)A = 
UL — 1104 oe ae on = U1 — riw x ù — EA, 
此 时 依旧 有 ft . 同样 的 , 需要 考虑 是 否 滑 动 ， 
有 : Fy = ee ae AN 考虑 力矩 , 有 : M 
= hô x Fr, HIN AL, ~ d* 
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2.4 引力 模型 a(t + dt) = aa (t+ ôt), č GE - a 
考虑 颗粒 之 间 的 引力 , 可 以 直接 在 颗粒 之 间 进 

行 计算 . 仍然 假设 颗粒 1 和 颗粒 2, 采用 和 上 面相 同 第 ee Midis V(t + 学) 和 角速度 矢量 


的 标记 , 定义 颗粒 1 和 颗粒 2 的 质量 分 别 为 mj 和 mo， 


有 颗粒 1 所 受到 的 颗粒 2 的 引力 : 

Fo = eae zı). 
考虑 计算 开销 , 假设 系统 由 "个 仿真 颗粒 组 成 , 那么 
每 次 计算 整个 系统 的 内 部 引力 作用 将 需要 计算 ml 
次 颗粒 之 间 的 两 两 引力 作用 , 计算 开销 巨大 , 因此 


我 们 需要 对 颗粒 按照 其 所 在 的 区 域 进行 分 块 计算 ， 
其 体 的 计算 细节 将 在 下 面 的 章节 中 展开 . 


2.5 ”积分 器 

在 数值 仿真 中 , 可 以 使 用 一 定 的 积分 器 数值 
算法 , 通过 已 有 的 动力 学 模型 六 = F(X,t) 和 to 时 
刻 初 值 , MAT AB AX O, 即 得 到 给 定 动力 
学 模型 和 初始 条 件 的 仿真 结果 . 出 于 对 于 精度 
和 计算 开销 的 综合 考虑 , 模型 里 采用 了 被 称 作 
“ 峙 跳 算 法 ”或 “Mid-Step Velocity Verlet” 的 积分 方 
[26] 78-82. 
第 一 步 : 将 位 置 矢量 元 
前 推算 一 步 (6b， 


im 


(t) 和 旋转 相位 矢量 Q(t) 


向 


+ U(t) dt + Taos ; 


O(t + dt) = O) + G(t) + Talar | 
其 中 , a(t) = “2, m 为 颗粒 的 质量 ; at) = ZO, 
7 为 颗粒 的 转动 惯量 . 
第 二 步 : 将 速度 矢量 Z(t) 和 角速度 矢量 w(t) 向 
前 推算 半 步 (区 )， 
(+ 号) = al), 
(t+ Ž) = alt) Laat 
第 三 步 : 计算 得 到 一 步 后 的 加 速度 a(t + 6t) 和 
角 加 速度 &(t + ôt), 


ak 


Ce er a(t + 60), a (t+ 


30-5 


如 此 便 可 以 以 积分 步 长 6t 将 仿真 向 前 推 


2.6 ”加 速算 法 
在 进行 多 颗粒 仿真 时 , 除了 考虑 模型 的 准确 1 
以 外 , 还 需要 考虑 计算 开销 的 问题 . 上 文中 我 们 提 
到 , 多 颗粒 系统 在 计算 时 有 极 大 的 计算 复杂 度 , 其 
中 , 对 计算 效率 影响 最 大 的 就 是 接触 算法 和 引力 算 
法 . 前 文 已 提 及 , 如 果 考 虑 颗粒 与 每 一 个 颗粒 之 间 
的 引力 相互 作用 , 那么 需要 进行 nl 次 颗粒 间 引 力 的 
计算 , 当 n 较 大 时 , 开销 巨大 , 几乎 无 法 进行 有 效 的 
计算 , 因此 , 仿真 过 程 需要 采取 一 定 的 近似 , 采用 合 
适 的 算法 来 对 仿真 进行 加 速 . 下 文 将 分 别 从 算法 、 
软件 和 硬件 层面 介绍 本 工作 采用 的 仿真 加 速 技 术 . 
2.6.1 “接触 加 速算 法 
在 进行 接触 力学 计算 时 , 需要 计算 两 颗粒 之 间 
ee eM A, KE > 0 时 , 认为 两 颗粒 发 生 碰撞 ， 


要 采用 上 文中 所 介绍 的 接触 力学 算法 对 颗粒 的 
受 力 进行 计 . 考虑 到 , d = |Z, 一 各 | 为 必须 进行 


计算 的 量 , 而 绝 大 多 数 颗 粒 都 不 与 某 一 给 定 的 颗 
粒 发 生 接 触 , 如 果 对 每 两 个 颗粒 都 进行 距离 的 计 
T, 则 产生 较 大 的 计算 开销 . 因为 某 一 给 定 的 颗粒 
仅 有 可 能 与 其 临近 的 颗粒 发 生 接 触 , 所 以 可 以 对 仿 
真 区 域 进行 网 格 划 分 . 不 失 一 般 性 地 考虑 矩形 六 面 
体 仿 真 区 域 , 可 将 其 划分 为 一 定 密度 的 正方 体 网 
格 区 域 , 仅 考 虑 一 个 颗粒 与 其 所 处 网 格 及 相 邻 网 
格 中 其 他 颗粒 的 接触 , 这 种 方法 被 称 为 Link List 方 
法 [261216 一 254. 

假设 仿真 区 域 为 zainyZmax] X [ymin Ymax] x 
[zmin; zamax] 所 围 成 的 矩形 六 面体 区 域 , 于 是 有 各 方 
向 上 的 长 度 : dj = Lmax — Lmin dy = Ymax — Ymin` 
ds 二 2max 一 Zmin, 正方 体 网 格 边 长 为 di, 由 此 得 到 
HAS Su de = 的 划分 好 的 方 格 区 域 . 为 了 满足 上 面 
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所 描述 的 某 个 颗粒 只 有 可 能 与 其 本 网 格 内 或 相 邻 
网 格 中 的 颗粒 发 生 接 触 的 条 件 , 有 四 = 2.5rmax, 其 
中 rmax 为 尺寸 最 大 的 颗粒 的 半径 . 

对 仿真 区 域 进 行 网 格 划分 之 后 , 可 以 对 每 一 个 
网 格 进行 标号 , M (is, ty, i.) RRA, 其 中 , 0 < is < 
d, 0 <S iy S 学 ,0 < iz < 他 ,由 此 可 以 对 所 有 网 
格 进行 遍历 , 并 进行 颗粒 间接 触 力 的 计算 . 注意 到 ， 
两 颗粒 之 间 的 接触 力 只 需要 计算 一 次 , 为 避免 重复 
计算 , 需要 采用 一 定 的 计算 策略 , 从 而 进一步 优化 
算法 . 


对 于 遍历 到 的 网 格 (i;, 纪 ,is), 如 图 3 所 示 , 从 上 
到 下 展示 了 上 、 中 、 下 3 层 的 网 格 分 层 结构 , 蓝 色 
方 格 为 网 格 (i,, 记 ,is), 绿色 为 此 时 参与 计算 的 网 
格 , 红色 为 不 参与 计算 的 网 格 , 以 此 遍历 所 有 的 网 
格 , 即 完 成 了 所 有 颗粒 的 接触 计算 . 由 于 在 接触 计 
算 时 得 到 了 颗粒 间 的 距离 , 此 时 同时 进行 颗粒 间 的 
引力 计算 将 节约 算 力 . 

此 算法 不 仅 可 以 保证 不 重复 遍历 网 格 , 即 只 通 
过 一 次 遍历 便 可 以 不 重复 地 计算 完 所 有 的 可 能 
生 的 颗粒 间接 触 , 另 一 方面 ,该 算法 也 可 以 保证 在 
并 行 计算 时 尽量 少 地 出 现 线程 竞争 的 情况 , 进一步 
提升 了 计算 效率 . 
2.6.2 ”引力 加 速算 法 

上 文中 提 到 , 如 果 考 虑 颗粒 与 每 一 个 颗粒 之 间 
的 引力 相互 作用 , 那么 需要 进行 nl 次 颗粒 间 引 力 的 
计算 , 当 n 较 大 时 , 开销 巨大 , 因此 , 需要 进行 一 定 
的 近似 和 算法 调整 来 优化 和 加 速 计算 过 程 271. 对 
于 某 一 颗粒 , 考虑 与 其 距离 较 远 的 颗粒 的 集合 , 它 
们 对 于 该 颗粒 的 引力 作用 可 以 近似 为 这 一 组 颗粒 
的 质心 处 一 质量 为 这 些 颗粒 总 质量 的 虚拟 颗粒 的 
引力 作用 . 在 这 种 假设 下 , 将 一 定数 量 网 格 的 集合 
区 域 组 织 为 正方 体 方 格 , 称 为 Cell. Cell 的 边 长 尺寸 
为 网 格 边 长 尺寸 的 整数 倍 , 并 充满 整个 仿真 区 域 
设 Cell 以 网 格 边 长 为 长 度 单位 的 边 长 尺寸 为 N., 一 
个 Cell 内 包含 的 网 格 数量 即 为 N3. 

对 于 同一 Cell 内 的 颗粒 , 在 计算 中 严格 计算 颗 
粒 间 的 引力 , 而 对 于 不 同 Cell 内 的 颗粒 之 间 的 引力 
作用 , 仪 考虑 当前 颗粒 与 其 他 Cell 内 处 在 质心 位 置 
的 质量 为 Cell 中 颗粒 总 质量 的 虚拟 颗粒 的 引力 作 
J. 在 计算 同一 Cell 内 的 颗粒 间 的 相互 引力 作用 时 ， 


全 


a 


如 图 4 所 示 为 在 一 个 Cell 内 进行 网 格 遍 历 的 情况 ， 
Ne = 6. 由 于 在 前 面 计算 接触 力 时 已 部 分 考虑 了 粒 
子 间 的 引力 ( 即 每 个 网 格 区 域 的 所 有 接触 网 格 都 已 
考虑 ), 因此 只 需 考 虑 该 Cell 内 尚未 考虑 的 晶 格 内 的 
粒子 即 可 . 蓝 色 网 格 为 此 时 遍历 到 的 网 格 , 红色 网 
格 表示 不 参与 计算 的 网 格 , 绿色 区 域 表示 参与 计算 
的 网 格 . 如 此 便 可 完成 Cel 内 所 有 颗粒 的 引力 计算 . 


3 ”当前 网 格 中 颗粒 与 临近 网 格 颗粒 碰撞 时 的 计算 策略 . 每 一 个 正方 
立方 块 代表 一 个 网 格 , 蓝 色 网 格 代表 当前 遍历 到 的 网 格 (is, iy, iz). 
能 与 该 网 格 中 的 颗粒 发 生 接触 的 颗粒 只 可 能 位 于 其 临近 网 格 中 , 即 
辐 所 示 的 上 下 左右 共 27 个 网 格 (包括 当前 网 格 本 身 ). 出 于 计算 效率 
E, 对 当前 遍历 到 的 网 格 (i;, iy, iis) 只 考虑 其 本 身 和 周围 26 个 网 格 
的 一 部 分 中 颗粒 可 能 发 生 的 接触 . 此 时 计算 考虑 的 网 格 用 绿 
色 标 出 , 不 参与 计算 的 网 格 用 红色 标 出 . 


TEFEN 


HH ii 


Fig.3 Calculation principle for collisions between grains in 
the current lattice and grains in its adjacent lattices. Every 
single cube in this figure represents a lattice, and the blue 
one represents the current lattice (iz, iy, 12) we calculate. 
Possible contacts only happen in adjacent lattices, and they 
are the 27 lattices (including the current lattice itself) 
plotted in this figure. Considering computational efficiency, 
only parts of the 26 adjacent lattices are involved in the 
calculation for the current lattice (is, tiy, iz). Lattices 
involved in the current calculation are marked with green 


while lattices not involved are marked with red. 


对 于 不 同 Cell 间 的 引力 , 由 于 需要 使 用 Cell 的 
质心 位 置 和 Cell 内 颗粒 总 质量 , 在 对 系统 进行 外 推 
时 每 一 步 均 计算 质心 位 置 和 总 质量 的 变化 并 存储 ， 
从 而 进一步 加 速 计算 
2.6.3 并行 算法 

随 着 计算 机 科学 技术 的 不 断 发 展 , 高 性 能 计算 
正成 为 解决 实际 问题 的 重要 工具 , 而 并 行 计算 作为 
高 性 能 计算 中 重要 的 细 分 领域 , 正在 计算 、 仿 真 、 
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数据 处 理 等 多 方面 发 挥 越 来 越 重要 的 作用 . 硬件 架 
构 上 来 看 , 并 行 计算 发 展 早期 主要 使 用 计算 机 中 央 
Ath FE AS (Central Processing Unit, CPU) 进 行 并行 ， 
于 CPU 设计 上 通用 性 的 需求 , 其 基本 计算 模块 的 
数量 较 少 , 图 形 处 理 器 (Graphics Processing Unit, 
GPU) 在 这 种 情况 下 诞生 , 成 为 并 行 计算 硬件 架构 
的 未 来 趋势 . 


图 4 计算 同一 Cel 内 颗粒 间 引 力 遍 历 网 格 时 考虑 的 网 格 示意 图 . 对 于 
司 一 Cell 内 的 颗粒 ,根据 算法 规定 需要 计算 颗粒 两 两 之 间 的 引力 作用 . 
出 于 计算 效率 的 考虑 , 在 Cell 内 遍历 网 格 时 , 对 于 当前 遍历 到 的 网 格 
HEE), 我 们 计算 引力 时 并 不 考虑 所 有 Cell 内 的 网 格 ,而 是 只 考虑 其 中 
的 一 部 分 . 图 中 考虑 的 网 格 用 绿色 标 出 , 不 参与 计算 的 网 格 用 红色 标 出 . 


Fig.4 Ca 


between grains in the current lattice and grains in its 


culation principle of lattice traversal for gravities 


neighboring lattices within the same Cell. For grains in the 
same Cell, gravities from each other are calculated based on 
the algorithm. Considering computational efficiency, when 
doing lattice traversal within a Cell, for the current lattice 
(marked with blue) we only calculate gravities from grains in 
parts of the lattices instead of all lattices. Lattices involved 
in the current calculation are marked with green while 


lattices not involved are marked with red. 


lhttps://www.openmp.org 
2https://www.open-mpi.org 
3https://developer.nvidia.com/cuda-toolkit 
fhttps://www.amdapp.org/s/ 
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从 并 行 计算 框架 上 来 看 , 目前 有 诸多 成 熟 


的 并 行 框架 可 供 选 择 . Open Multi-Processing 
(openMP)', Open Message Passing Interface 
(openMPI)? 为 代表 性 的 适用 于 不 同 计算 环境 和 
条 件 的 CPU 并 行 计算 框架 ,C+ 十 11 之 后 也 提供 了 
<thread> 标 准 库 供 并 行 计算 使 用 , 使 用 CPU 并 行 
能 简单 快速 地 将 程序 部 署 到 服务 器 上 ， 从 而 利 
用 服务 器 的 多 核 计 算 资 源 加 快 计 算 速 度 . 此 外 ， 
NVDIA 和 AMD 公 司 也 相应 开发 了 基于 自身 生产 
的 GPU 的 并 行 计算 框架 Compute Unified Device 
Architecture (CUDA)? Accelerated Parallel Pro- 
cessing (APP)*, 通过 使 用 个 人 电脑 的 图 形 处 理 器 ， 
可 以 方便 地 将 并 行 计算 程序 部 署 在 本 地 . 
本 工作 中 , 使 用 了 C++ 中 <thread> 标 准 库 提 
供 的 并 行 计算 框架 解决 方案 , 目的 在 于 利用 C 十 十 和 
其 标准 库 便于 移植 的 特性 , 开发 可 移植 、 容 易 部 署 
的 并 行 计算 程序 . 程序 中 不 同 的 力 模 型 的 计算 
由 CPU 串 行 调用 不 同 的 力 模 型 计算 函数 , 函数 内 部 
采用 并 行 结构 , 以 Cell 为 基本 并 行 单元 , 灵活 地 变 
更 并 行 线程 的 结构 , 从 而 增加 计算 效率 . 同时 , 根据 
前 文 所 述 , 采用 了 合理 的 Cel] 或 网 格 遍历 策 略 , 使 
线程 竞争 最 小 化 , 进一步 提升 并 行 效率 . 
如 图 5 所 示 , 分 别 用 点 标示 了 单线 程 和 多 线程 
(4 线程 ) 测 试 的 结果 , 测试 了 不 同 数量 测试 颗粒 的 
数值 实验 所 耗费 的 时 间 . 测试 总 时 间 为 tiot, 测试 颗 
粒 总 数 为 n. 使 用 函数 tiot = ksn 对 单线 程 结果 进行 
了 拟 合 , 绘制 了 红色 直线 , 随后 绘制 直线 ttot = 笠 


4 


与 多 线程 的 预期 结果 符合 得 很 好 , 并 行 效率 很 高 . 
2.7 ”仿真 参数 的 选取 

在 仿真 中 , 除了 选择 合适 的 积分 器 来 减 小 仿真 
误差 之 外 , 更 重要 的 是 对 于 积分 器 步 长 参数 和 仿真 
模型 参数 的 合理 选取 . 在 本 工作 中 , 主要 讨论 法 向 
弹性 系数 血 、 阻 尼 系 数 Y 和 仿真 步 长 %t 的 选取 , 在 
仿真 过 程 中 , 6: 和 和 操 共 同 决定 了 仿真 的 外 推 精度 ， 
7 决定 了 能 量 耗 散 的 速度 . 


n, 
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2.7.1 如 的 选取 

一 方面 , 忽略 耗 散 系数 , 即 y = 0, 考虑 两 个 正 
向 碰撞 的 小 球 ， 有 能 量 守恒 : imav = 3kné2ax， 
其 中 mef = OM, vo 为 接触 前 两 个 颗粒 之 间 的 相 
对 速度 ， twax 为 接触 过 程 中 颗粒 的 最 大 “形变 ” 
在 软 硬 球 模型 中 , 我 们 希望 有 较 小 的 “形变 ”, 于 
fe i MRK NAMA, 根据 具体 的 仿真 条 件 ， 可 
DIN BAe Wk, fA. 在 非 线 性 弹 千 模 型 中 , 非 线 性 


[fini 
uy 


学 报 3A 
E= 了 exp (-;-t) E cos (wt) + > sin (wt)] ; 
Eh, w= /起 -一 Ta 于 是 整个 碰撞 发 生 的 时 
间 为 : 

(223 7 (1) 
O a kn y? 
Mort 4 Me 


WL) REAR, “ky BERING AY ERRE, 
此 时 若 要 求 积分 过 程 中 有 较 多 的 积分 次 数 , 积分 步 
长 需要 取得 很 小 , 仿真 的 总 体 计算 开销 便 会 增 大 


弹性 系数 kr 中 值 可 由 实验 测量 的 材料 泊 松 比 o、 
杨 氏 模 量 已 和 接触 颗粒 的 半径 r 来 计算 叶 ， 有 : 


(non) _ 4_E 
hypo) = Sav. 
5000 T T 
e 4-thread 
e 1-thread 
4000 F . 


3000r 


ttot / S 


2000} 


1000 F 


000 1250 1500 1750 
n 

色 和 蓝 色 的 数据 点 分 别 标示 了 单线 
函数 tot = ksm 对 单线 程 结果 


学 n 的 


2000 


15 并 行 计算 效率 测试 结果 . 用 红 
程 和 多 线程 (4 线程 ) 测 试 的 结果 . 使 
行 了 拟 合 , 得 到 了 红色 拟 合 直 线 , 同时 根据 拟 合 结果 画 出 ttot 
直线 图 , 与 多 线程 的 结果 预期 符合 得 很 好 , 并 行 效率 很 高 . 


Fig.5 The test result of the parallel computation efficiency. 
1-thread and multi-thread (4-thread) are marked 


Results o 
with red and blue scatter points respectively. Function 
ttot = ksn is used for a fitting of the 1-thread result and the 
Ea n is 


red line is obtained. A blue line representing ttot 


plotted based on the 4-thread result, which meets the 


multi-thread computation result well and the parallel 


computation efficiency is high. 
另 一 方面 , 在 接触 的 积分 过 程 中 , 我 们 希望 有 
较 多 的 积分 次 数 , 而 铝 值 决定 了 接触 的 持续 时 间 . 
容易 写 出 的 常 微分 方程 : 
mené T y +k = 0. 
这 是 一 个 二 阶 常 系数 微分 方程 , 由 初 值 : t = 0, 
E = 0; = wo 可 以 给 出 解析 解 : 


因此 , 需要 选择 合适 的 总 值 , 从 而 达到 模型 合理 和 
计算 效率 的 平衡 . 

Schwartz 等 上 癌 在 描述 DEM 仿 真 程序 PKDGRAV 
时 考虑 了 对 总 选取 的 两 种 情况 . 第 一 种 情况 考虑 动 
态 的 颗粒 碰撞 接触 过 程 , 在 采用 线形 弹簧 模型 的 前 
提 下 , 可 以 由 wm 及 Snax 得 到 估计 的 总 值 . 由 能 量 守 
ER mavi = 3kn62。x 容 易 得 到 ks ~ mea EB )?. 
第 二 种 情况 考虑 低速 的 颗粒 聚集 相互 挤 压 形成 的 
准 静 止 状态 , 通过 佑 计 颗 粒 再 接触 过 程 中 的 最 大 转 
FRA Bll ey (EL. 假设 高 度 瓦 的 容器 中 装 满 半 径 为 s 的 
小 球 , 小 球 密度 为 p, 在 重力 加 速度 ao 作用 下 底层 
小 球 受 到 ~ gpagHHs? 的 压力 , 其 中 9 为 小 球 颗粒 
WARE, 于 是 pp 为 系统 总 体 密度 , 因此 有 ,~ 
teas” ,在 实际 仿真 过 程 中 选取 记 值 可 以 参考 二 
者 之 间 的 较 大 值 . 
2.7.2 7 的 选取 


将 (1) 式 得 到 的 时 间 代 入 方程 的 解 , 可 以 得 到 
反弹 之 后 两 微粒 之 间 的 相对 速度 : 
TY 
v= wo ( Vi me a) 
速度 的 衰减 由 耗 散 系数 y 决 定 , 于 是 有 : 
voo TY 
eed 
这 一 数值 可 以 由 实验 测量 得 到 . 
2.7.3 6t 的 选取 
6t 的 选取 决定 了 仿真 外 推 过 程 中 积分 器 带 来 
的 误差 的 大 小 , 在 颗粒 接触 过 程 中 , 积分 的 步 数 
越 多 , 结果 越 准 确 . 一 个 直观 的 想法 是 使 用 上 面 的 


Ty 


gE ) 来 进行 5 的 选取 . 考 


指标 地 = exp (- 
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处 (1) 式 给 出 的 接触 持续 时 间 , 将 积分 步 长 6 表示 为 

如 下 形式 : 

f=C =C 
W 


T 


= 0) 


Meff 4m > ff 


其 中 C 为 一 个 常数 . 在 仿真 之 前 , 可 以 根据 仿真 条 
件 , 计算 一 个 简单 算 例 , 考虑 仿真 中 可 能 出 现 的 vo 
最 大 的 情形 下 , 当前 C 值 所 决定 的 世 值 是 否 与 理论 
分 析 相 当 , 如 果 相 差 较 大 , 则 对 C 进 行 更 改 并 迭代 ， 
直到 蕊 符合 理论 分 析 的 值 . 

事实 上 , 在 绝 大 部 分 情况 下 , StH PAE H ky 
定 , 因为 阻尼 系数 7 在 大 多 数 情况 下 较 小 , 则 有 (2) 


式 的 近似 : öt ~ Cr TE. 数值 实验 表明 , 对 于 相 
同 的 如 值 , 不 同 的 7 取 值 时 , 即使 始终 采用 y = 0 时 
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对 于 给 定 的 各 和 一 组 y 值 , 使 用 7 = 0 人 迭代 获得 C 值 ， 
对 所 有 相同 加 下 的 y 采 取 同 一 C 值 , 即 同一 积分 步 
长 进行 计算 . 取 两 个 颗粒 密度 p = 2 x 103 kg.m-3、 
颗粒 半径 x = 0.1 m 的 球形 颗粒 单元 , 发 生 相 向 初 
始 相 对 速度 wo = 0.01 m- s-! 的 碰撞 , 记录 碰撞 
后 的 相对 速度 v, 然后 计算 仿真 结果 (esm = 之 ) 的 
值 , 于 是 可 以 得 到 上 述 颗 粒 条 件 下 的 各 、7 和 仿真 
耗 散 es = 芝 的 关系 图 . 结果 如 图 7 所 示 , 其 中 白 
色 部 分 表示 4megks — 7? < 0, 为 过 阻尼 情况 ， 
进行 仿真 , 其 余部 分 颜色 与 颜色 条 上 不 同 各 取 值 
相对 应 , 表征 当前 如、?7 取 值 下 仿真 得 到 的 系统 
能 量 耗 散 结果 ( 左 图 ), 与 相应 的 理论 值 做 差 , 计算 
|esim 一 etheory| 得 到 理论 值 和 仿真 值 的 差 值 ( 右 图 )， 
TY 
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的 C 值 进行 计算 , 所 得 到 的 结果 均 能 保证 基本 与 二 
体 接触 问题 中 的 世 理 论 值 相 符 . 


3 测试 算 例 
3.1 已 、yY、6t 的 取 值 实验 

EXTA Tk. y StH RUB, 下 面 通 过 
实验 来 验证 上 文 的 理论 推导 和 分 析 . Ak, = 10° N. 
m-1, y = 100 N.s.m-1, 为 不 失 一 般 性 , 考虑 一 颗 
粒 与 平面 的 接触 , 颗粒 从 一 定 高 度 下 落 , 经 平面 反 
弹 后 达到 新 的 高 度 , 由 于 耗 散 的 存在 , 新 的 高 度 较 
原 高 度 较 低 . 取 颗 粒 密 度 p = 2 x 103 kg . m-?,， 颗 
粒 半径 均 为 r = 0.1 m, 此 时 有 e = 2 ~ 0.84, 容易 
推导 , 前 后 的 高 度 ho 与 的 关系 为 : + = VE, 取 初 
AC 二 0.1, 如 果 仿 真得 到 的 R/E 与 理论 值 0.84 
相差 过 大 , 则 认为 此 时 的 积分 步 长 选取 不 合适 , C 
> S, 如 此 和 迭代, 得 到 合适 的 积分 步 长 丝 , 如 图 6 所 
示 ， 选择 之 或 V 声 作为 参考 指标 ， 当 C 值 足够 小 时 ， 
二 者 均 能 够 收敛 到 理论 值 . 
进一步 , 我 们 可 以 计算 不 同和 > 值 下 的 耗 散 
车 的 理论 值 , 同时 针对 不 同 的 各 和 7 值 仿真 得 到 的 
结果 并 与 理论 值 进行 比较 . 上 文中 提 到 , kate t fih 
过 程 中 起 到 决定 性 作用 , 因此 对 于 相同 的 如 在 不 
同 的 y 取 值 下 , 可 以 使 用 相同 的 C 值 , 即 相 同 的 积 4 
步 长 . 图 7 所 示 的 数值 实验 即 采取 这 样 的 计算 策略 ， 


TS 


其 中 ，etheory = exp ( as): 由 图 可 知 仿 


真 结果 esim 与 前 述 理论 的 etneory 符 合 得 较 好 . 


— eH” 


而 
二 ,证 
= e=Vt 


107? 1071 


È 


图 


16 迭代 C 值 得 到 合适 的 积分 步 长 . 图 中 
化 e 值 的 变化 情况 . 数值 试验 的 结果 用 蓝 
e= thle = /总 两 种 方式 仿真 得 


线 代 表 随 着 C 值 的 变 
色 和 检 色 分 别 表示 使 
e 值 . 随 着 C 值 的 减 小 , 两 种 


数值 仿真 结果 均 能 够 收敛 到 理论 的 e exp ( 


T ). 
Fig.6 C value iteration for suitable time step for 
propagation. The curves plotted are e values with respect to 
the change of C value. Numerical results are plotted in blue 


and orange for simulated e values using function e = aa and 


function e = 


e = exp ( 


+ respectively. The two numerical results 
0 
are able to converge to theoretical value 


TY 
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) with the decrease of the C value. 
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esim 一 etheory| 
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y/(N-m-s7) 


Ka 


v 
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图 


7 通过 数值 实验 仿真 得 到 的 k&,。、Yy、 世 关系 


(£ 


度 vo = 0.01 m . s 的 磁 撞 后 的 相对 速度 v 得 到 , 定义 有 仿真 的 esim = 2. 数值 仿真 的 结果 与 理论 值 的 差 值 


ny 


和 数值 仿真 结果 与 理论 值 的 差 什 


图 


(E). 数值 试验 的 结果 通过 测量 两 颗粒 发 生 初 始 相 对 速 


H etheory exp ( 


ae 
vo 


Fig.7 Simulation results of kn, y, 
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relationship contour (left) and the deviation of simulation value and theoretical value 


). 


计 |€sim etheory| 得 到 ， 其 


仿真 结果 与 理论 值 符合 得 很 好 . 


(right). The results of numerical tests come from measuring the relative velocity v after a collision happens between two grains 


with an initial relative velocity of vp = 0.01 m-s~' 


simulation values and theoretical values are calculated using |esim 一 etheory| and etheory exp ( 


and the simulated result is defined as esim = Bo" The deviation between the 


TY 


V4mergkn =y? 


) . The simulation 


result meets the theory well. 


这 一 实验 结果 说 明了 选取 的 迭代 算法 能 够 有 
效 地 在 程序 中 确定 适当 的 步 长 , 且 随 着 积分 步 长 的 
缩短 , 得 到 的 结果 应 更 准确 . 进而 通过 对 恕 、7 和 耗 
散 艺 关系 的 仿真 说 明了 此 多 颗粒 仿真 算法 及 程序 
在 两 颗粒 间 及 颗粒 与 平面 接触 的 问题 上 能 够 得 到 
合理 、 准 确 的 结果 , 有 效 性 得 到 了 验证 . 


3.2 ”声波 速度 的 仿真 测量 

THESE RS, 声波 在 碎 石 堆 中 的 传播 速度 
是 重要 的 可 观测 参数 , P 波 和 8 波 在 碎 石 堆 中 的 传 
播 速度 与 组 成 颗粒 单元 的 物理 性 质 紧 密 关 联 , 因 
此 , 对 于 声波 传播 速度 的 研究 对 宏观 仿真 体 的 内 部 
结构 研究 、 物 理性 质 研 究 有 重要 的 意义 . 实验 研究 
KAHOL, 硅 酸 盐 或 玻璃 材料 碎 石 堆 中 的 声速 约 为 
170 m- sot, 我 们 可 以 使 用 本 软件 来 仿真 这 一 声速 
数据 . 

设置 12.5 mm x 12.5 mm x 50 mm 的 仿真 区 域 ， 
底部 和 顶部 采用 接触 边界 , 即 对 应 于 > = 0 mm 和 


mh 
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z = 50 mm 设置 接触 边界 条 件 , 四 壁 采用 循环 边 
界 条 件 进 行 构造 . 在 仿真 区 域 中 添加 颗粒 (直径 
为 0.5 mm 的 均匀 球 ), 基于 颗粒 尺寸 选择 边 长 为 
0.625 mm 的 网 格 尺寸 , 对 应 此 例 中 区 域 被 划分 为 20 
x 20 x 80 的 网 格 区 域 , 选取 Cell 的 边 长 尺寸 为 5 个 
网 格 边 长 . 

选择 颗粒 的 密度 p = 2500 kg - m-3 为 一 般 硅 酸 
盐 的 密度 . 从 仿真 区 域 底部 开始 , 放置 所 需要 的 实 
验 颗 粒 , 一 直到 仿真 区 域 的 颗粒 高 度 达 到 18 mm 为 
ik. 

使 用 的 力学 模型 包括 接触 力 和 用 于 模拟 重力 
环境 的 定常 力 . 采用 两 组 接触 系数 进行 计算 , 接触 
系数 的 选择 参考 了 Sinchez 等 2 相同 数值 实验 中 的 
取 值 , 便于 对 实验 结果 进行 比较 , 如 表 1 所 示 . 使 用 
Leap Frog 单 步 法 积分 器 进行 积分 , 积分 格式 向 前 
积分 , 初始 设置 步 长 为 5 x 10-9 s, 步 长 随 着 力 模型 
参数 的 选择 通过 公式 5 = Cf met 自主 选择 . 
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表 1 2 组 声波 传播 实验 参数 选取 表 


Table 1 Parameters table of two sound wave propagation numerical experiments 


Case Index kn/(N-m7') k/(N-m7') y/(N-s-m7") p 
1 10° 2.8 x 104 71x 107% 0.5 
2 1.8 x 10+ 5.1 x 10° 3x 107% 0.5 


首先 对 整个 仿真 系统 进行 静 置 稳定 , 在 不 同 
的 重力 加 速度 值 下 , 系统 的 能 量 逐 渐 趋 于 定 值 ， 
此 时 对 位 于 底部 的 颗粒 施加 一 个 z 方 向 大 小 为 0.5 
m.s-! 的 速度 增 量 , 这 一 速度 增 量 的 影响 会 在 整 
个 系统 中 传播 . 我 们 可 以 通过 仿真 这 一 增 量 传 择 
的 速度 来 获得 这 一 系统 中 的 声波 的 速度 . 如 图 8 所 
示 , 两 组 不 同 颗粒 参数 的 实验 结果 分 别 使 用 三 角 和 
形 散 点 表示 . 在 实验 中 对 同一 组 颗粒 参数 , 取 不 
同 的 重力 加 速度 g 值 (以 地 球 表 面 重力 加 速度 g。 = 
一 9.8 m :s 悦 为 单位 ), 测算 出 不 同 深度 处 的 声波 
播 速度 . 两 组 数值 实验 的 结果 都 表现 出 类 似 的 
A, 随 着 声波 的 传播 深度 的 增加 , 声速 逐渐 趋 于 
个 定 值 , 这 一 值 与 重力 加 速度 g 值 的 关联 并 不 明显 ， 
但 显著 地 受到 接触 系数 的 影响 . 第 二 组 实验 所 选取 
的 接触 系数 导致 仿真 系统 中 的 声速 逐渐 趋 近 于 硅 
酸 盐 或 玻璃 材料 碎 石 堆 中 的 声速 , 即 红 线 标示 出 
的 170 ms-1. 


3.3 RAE Te 

对 于 碎 石 堆 结构 的 小 行星 , 可 以 使 用 离散 元 方 
法 对 其 进行 仿真 . 选取 1000 个 半径 为 r = 1 m, 密度 
=3x103kg.m-3 的 仿真 颗粒 单元 , 将 仿真 颗 
粒 随机 分 布 在 50 m x 50m x 50 m 的 仿真 区 域 中 ， 
图 9 所 示 为 颗粒 的 初始 状态 分 布 情况 . 设置 力学 模 
型 , 考虑 引力 和 接触 力 的 作用 , 采用 接触 系数 : 法 向 
弹性 系数 户 = 1.8 10°N-m!, WAAR, 
= 2k,, Mitt Rž = 3 x 10'N-s-m!, BRA 
数 / = 0.5. 接触 系数 的 选择 参考 了 S&nchez 等 Bol 针 
对 较 大 颗粒 尺寸 的 小 行星 设置 的 接触 系数 . 对 于 较 
小 的 颗粒 尺寸 和 密度 , 使 用 这 一 组 接触 系数 使 得 稳 
定 小 行星 中 每 一 组 接触 颗粒 有 更 小 的 嵌入 深度 . 对 
系统 进行 静 置 , 仿真 颗粒 逐渐 聚集 成 团 , 20000 s 后 
基本 达到 静止 , 其 状态 如 图 10 所 示 , 考虑 初始 颗粒 
分 布 的 状态 为 近似 均匀 分 布 在 一 个 矩形 多 面体 空 


E. 
E 


E 
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间 内 , 最 终 聚 合 为 一 个 类 似 多 面体 的 结构 . 


2000 F (Case 1) 9g9= 10-29e。 | 
(Case 1) g= 10-1ge 
1750F (Case 1) g=10°g. | 
(Case 1) g= 1019。 
T 1500 (Case 1) g=10?ge | 
E eer (Case 2) g=1072g. | 
三 (Case 2) g= 10-1ge 
D =10° 
® 1000 (Case 2) g= 10°ge | 
a (Case 2) g=101g, 
3 750} (Case 2) g=10?ge | 
wn 
500 
2507 
0.0000 0.0025 0.0050 0.0075 0.0100 0.0125 0.0150 0.0175 
Height /m 
图 8 声速 仿真 结果 图 . 两 组 不 同 颗粒 参数 的 实验 结果 分 别 使 用 三 角 和 
到 形 散 点 表示 . 在 实验 中 对 同一 组 颗粒 参数 , 取 不 同 的 重力 加 速 
度 g 值 (以 地 球 表面 重力 加 速度 g。 = 一 9.8 m-s 7 AMAL), 测算 出 不 


同 深度 处 的 声波 传播 速度 . 
Fig.8 Results of the sound wave propagation numerical 
experiment. Results of two sets of particle parameters are 
scattered in triangles and circles. For each set of particle 
parameters, sound wave propagation velocity is measured at 
different height using different gravitational acceleration g 
values (with the gravitational acceleration on the Earth’s 


surface ge = —9.8 m - s7? as a unit). 


fe POR, 以 距离 碎 石 堆 质心 的 距离 为 标准 , 保 
留 距离 4 < 10 m 的 仿真 颗粒 , 继续 静 置 仿真 系统 ， 
在 这 一 过 程 中 , 系统 逐渐 趋 于 稳定 , 得 到 近似 球体 
的 、 半 径 约 10 m 的 小 行星 , 如 图 11 所 示 . 

对 密度 均匀 的 球形 小 天 体 , 其 内 部 一 点 的 压强 
计算 公式 为 P = enGp, (he —r7), HARA) 
的 半径 , ”为 小 行星 质心 到 小 行星 内 某 一 点 的 吕 
而 对 于 我 们 的 仿真 系统 , 有 


Ni No 


1 oe 
Pe = 5 2 Sifi 


i=0 j=0 
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FEAR AER IG ER TE AY ba, 不 同 标号 标示 选择 距离 
Wi» ASAE VERSE, 5; 为 第 k 个 球形 球 壳 的 表面 
BA, 和 Ni 为 该 球 壳 外 部 与 球 壳 k 接 触 的 颗粒 总 数 ，N2 
为 与 颗粒 ;接触 且 位 于 球 壳 内 部 的 颗粒 总 数 ， 记 yj 为 
颗粒 ;受到 颗粒 7 的 接触 力 , 多 为 颗粒 ?位置 指向 碎 石 
全 质心 位 置 的 单位 矢量 . 考虑 球形 颗粒 组 成 的 碎 石 
EWER, 有 小 行星 密度 pa = Zp, 图 12 给 出 了 理论 


值 与 仿真 值 的 比较 , 不 难看 出 二 者 符合 得 很 好 . 


对 


9 初始 时 刻 的 仿真 颗粒 分 布 . 1000 个 半径 为 r = 1 m 的 颗粒 随机 分 
布 在 50 m x 50m x 50 m 的 仿真 区 域 中 . 


Fig.9 The original grain distribution. 1000 grains with 
radius r = 1 m are randomly distributed in a 50m x 50m x 


50 m simulation area. 


图 


10 ”20000 s 后 的 仿真 颗粒 分 布 


Fig.10 Grain distribution after 20000 s simulation 


3.4 ”旋转 稳定 性 

小 行星 的 自转 是 影响 小 行星 演化 进程 的 重 
要 因素 ,，Yarkovsky-O’Keefe-Radzievskii-Paddack 
(YORP) 效 应 叫 、 大 天 体 的 近 距 离 飞 越 等 都 能 够 
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3 期 


影响 小 行星 的 自转 速度 . 对 双星 系统 而 言 , 两 星 
间 的 轨 旋 耦合 效应 、 潮 汐 效应 以 及 可 能 存在 的 Bi- 
nary YORP (BYORP) 效 应 都 会 改变 小 天 体 的 转 
速 , 目前 认为 YORP 效 应 是 造成 近 地 双 星 和 小 尺寸 
主 带 双星 的 主要 成 因 [B83. 本 工作 仅 考 虑 单 星 匀 速 自 
转 的 情形 , 在 不 考虑 粘 聚 力 的 前 提 下 , 当 有 自转 速 
度 超过 we = CM int, 小 行星 无 法 继续 保持 引力 和 
内 部 接触 压力 的 平衡 , 小 行星 开始 裂解 . 


图 11 静 置 后 得 到 的 近 球 形 小 行星 
Fig.11 Settled spherical asteroid 
0.035 
e 
e Simulation 
Theory 
0.030 
0.025 
æ 0.020 
a 
ĉ 0,015 
0.010 
0.005 
0.000, 4 6 
Distance From the Center/ m 
图 12 ”小 行星 内 部 压强 仿真 


Fig.12 Asteroid inner pressure scatter comparison between 


simulation and theory 


继续 使 用 上 一 个 数值 实验 中 所 计算 的 球形 小 
行星 模型 , 使 其 从 静止 开始 以 (V2, V2, 0) 为 轴 自 转 . 
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从 0.1 wc 的 角速度 开始 计算 , 以 0.1 w. 为 角速度 增 。 程 仿真 的 可 靠 性 , 证 明 程序 仿真 的 基本 颗粒 接触 过 


加 的 步 长 直到 w = we, 每 一 个 旋转 速度 计算 50 s. 程 是 可 靠 的 ; 

每 1 s 获 取 一 次 z2 十 好 = 102, 20 < z < 30 的 环形 区 第 3.2 节 中 声波 速度 的 仿真 测量 证 明了 碎 石 堆 
域 ( 即 赤 道上 ) 的 平均 每 个 颗粒 所 受 的 向 小 行星 外 部 环境 下 多 颗粒 相互 碰撞 带 来 的 接触 力学 复杂 耦合 
的 接触 压力 , 对 同一 个 旋转 速度 仿真 时 间 内 所 获取 计算 的 准确 性 和 有 效 性 , 是 对 后 续 碎 石 堆 小 行星 仿 
的 所 有 压力 值 取 平均 , 当 w = we 时 , 小 行星 此 时 环 真 可 靠 性 的 有 力 验 证 ; 

区 域内 颗粒 不 发 生 接 触 , 小 行星 发 生 破碎 裂解 . 第 3.3 节 中 碎 石 堆 小 行星 从 分 散 颗 粒 的 引力 凝 
如 图 13 所 示 , 图 中 展示 了 在 小 行星 旋转 的 不 断 加 速 聚 过 程 验证 了 程序 中 对 引力 的 近似 处 理 的 可 靠 性 


过 程 中 小 行星 赤道 上 环形 区 域 颗粒 平均 受到 接触 
压力 的 变化 情况 . 小 行星 在 逐渐 加 速 旋转 的 过 程 
中 , 接触 力 平衡 了 小 行星 的 引力 , 在 低速 转动 时 保 


持 了 
星 的 


接触 大 大 减少 或 不 接触 , 即 压 力 趋 于 0, 小 行星 发 生 


分 裂 . 
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的 构 
接触 
速 计 
软件 


整个 系统 的 平衡 , 当 旋转 速度 超过 we 时 , 小 行 
引力 无 法 维持 旋转 向 心力 , 此 时 组 成 颗粒 之 间 


和 


时 


有 效 性 , 其 内 部 分 层 压 强 的 数值 实验 测量 与 理论 


= 


值 
程 


了 


吻合 说 明 程 序 的 力学 模型 仿真 和 动力 学 外 推 过 
是 真实 物理 条 件 的 有 效 、 可 靠 的 仿真 ; 
第 3.4 节 中 碎 石 堆 小 行星 旋转 分 解 的 过 程 验证 
程序 仿真 小 行星 旋转 稳定 性 过 程 的 可 靠 性 , 为 后 


续 可 能 进行 的 数值 实验 提供 了 旋转 力学 的 验证 和 


参考 . 现 有 的 使 


j 离 散 元 方法 的 研究 中 , 潮汐 效应 、 


YORP 或 BYORP 效 应 对 于 小 行星 自转 的 研究 依赖 


— Theory 
e Simulation 


0 02 0.4 0.6 08 1.0 12 
Angular Velocity / we 


转速 度 与 内 部 环形 区 域 平均 接触 压力 的 关系 


13 ”小 行星 


13 Relationship between average force in an annular 


area and the asteroid’s spin speed 


总 结 


10^ 2A 


文章 首次 介绍 了 作者 开发 的 离散 元 仿真 软件 
建 细 节 , 包括 基本 仿真 单元 的 选取 、 粒 子 间 的 
力 建 模 和 引力 建 模 以 及 软件 中 采取 的 一 些 加 
算 措施 ; 其 次 , 通过 几 个 数值 仿真 实验 验证 了 
的 仿真 结果 的 可 信 性 , 这 些 实验 包括 : 


第 3.1 节 中 高 、y、6t 的 取 值 实验 验证 了 接触 过 
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m, 
而 


ab 
He 


验证 一 些 猜 测 也 是 作者 开发 此 软件 的 初衷 . 在 后 期 
的 仿真 模拟 中 , 作者 会 进 一 


于 基于 旋转 角速度 矢量 演化 预测 的 离散 元 仿真 过 

这 一 过 程 常 常 因 碎 石 堆 小 行星 形状 变化 的 影响 
与 刚体 力学 表现 出 显著 不 同 , 是 未 来 相关 研究 可 
的 重要 方向 . 对 这 些 过 程 进 行 仿真 模拟 以 提取 、 


步 完 善 现 有 的 软件 , 加 


入 不 同 基本 单元 的 选取 、 粘 聚 力 的 模拟 等 
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Discrete Element Method Simulation System for Asteroid 


LIU Mu-lin 


HOU Xi-yun 


(Department of Astronomy and Space Science, Nanjing University, Nanjing 210023) 


ABSTRACT Asteroid detection is now a hot spot of solar system exploration. The understanding of the 
evolution of asteroids has a great benefit on researches to the origin of the solar system. An important 
topic of the evolution research is the evolution of inner structures of asteroids, in other words, asteroids’ 
evolution of shapes and structures under multiple mechanisms. A common method to simulate the dy- 
namical evolution of asteroids is Discrete Element Method (DEM) under the assumption that asteroids 
are in rubble-pile structures. Some teams have developed several kinds of softwares for DEM simulations. 
The basics, realization, algorithms of our software “Multi-particle system simulation software based on 
DEM algorithm” are introduced in this article and the software is verified using two body collision, sound 
wave propagation, inner stress of asteroid, and asteroid spin disruption setups. 


Key words methods: numerical, asteroids: general, asteroids: evolution 
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